Pair contact process with diffusion — 
A new type of nonequilibrium critical behavior? 
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In the preceding article Carlon et al. investigate the critical behavior of the pair contact process with 
diffusion. Using density matrix renormalization group methods, they estimate the critical exponents, 
raising the possibility that the transition might belong to the same universality class as branching 
annihilating random walks with even numbers of offspring. This is surprising since the model does 
not have an explicit parity-conserving symmetry. In order to understand this contradiction, we 
estimate the critical exponents by Monte Carlo simulations. The results suggest that the transition 
might belong to a different universality class that has not been investigated before. 



Symmetries and conservation laws are known to play 
an important role in the theory of nonequilibrium critical 
phenomena Q. As in equilibrium statistical mechanics, 
most phase transitions far from equilibrium are charac- 
terized by certain universal properties. The number of 
possible universality classes, especially in 1+1 dimen- 
sions, is believed to be finite. Typically each of these 
universality classes is associated with certain symmetry 
properties. 

One of the most prominent universality classes of 
nonequilibrium phase transitions is directed percolation 
(DP) According to a conjecture by Janssen and 

Grassberger, any phase transition from a fluctuating 
phase into a single absorbing state in a homogeneous sys- 
tem with short-range interactions should belong to the 
DP universality class, provided that there are no special 
attributes such as quenched disorder, additional conser- 
vation laws, or unconventional symmetries QJ^]. Roughly 
speaking, the DP class covers all models following the 
reaction-diffusion scheme A «-> 2A, A — > 0. Regarding 
systems with a single absorbing state the DP conjec- 
ture is well established nowadays. However, even various 
systems with infinitely many absorbing states have been 
found to belong to the DP class as well [|| |) . 

Exceptions from DP are usually observed if one of 
the conditions listed in the DP conjecture is violated. 
This happens, for instance, in models with additional 
symmetries. An important example is the so-called 
parity-conserving (PC) universality class, which is repre- 
sented most prominently by branching annihilating ran- 
dom walks with two offspring A —> 3 A, 2 A —> [| 0. 
In 1+1 dimensions this process can be interpreted as 
as a ^-symmetric spreading process with branching- 
annihilating kinks between oppositely oriented absorb- 
ing domains. Examples include certain kinetic Ising 
models fl2|| , interacting monomer-dimer models [fl3|| , as 
well as generalized versions of the Domany-Kinzel model 
and the contact process with two symmetric absorbing 
states 1T1. 



A very interesting model, which is studied in the 
present work, is the (l+l)-dimensional pair contact pro- 
cess (PCP) 2A -> 3A, 2A -» Q. Depending on 
the rate for offspring production, this model displays a 
nonequilibrium transition from an active into an inactive 
phase. Without diffusion the PCP has infinitely many 
absorbing states and the transition is found to belong 
to the universality class of DP. The pair contact process 
with diffusion (PCPD), however, is characterized by a 
different type of critical behavior. In the inactive phase, 
for example, the order parameter no longer decays expo- 
nentially, instead it is governed by an annihilating ran- 
dom walk with an algebraic decay. Moreover, the PCPD 
has only two absorbing states, namely, the empty lattice 
and the state with a single diffusing particle. For these 
reasons the transition is expected to cross over to a dif- 
ferent universality class. The PCPD, also called the an- 
nihilation/fission process, was first proposed by Howard 
and Tauber jl6| as a model interpolating between "real" 
and "imaginary" noise. Based on a field-theoretic renor- 
malization group study, they predicted non-DP critical 
behavior at the transition. 

In the preceding article, Carlon, Hcnkcl, and 
Schollwock [[t7J investigate a lattice model of the PCPD 
with random-sequential updates. In contrast to Ref. |l6| ], 
each site of the lattice can be occupied by at most one 
particle, leading to a well-defined particle density in the 
active phase. Performing a careful density matrix renor- 
malization group (DMRG) study jl8|,[l9|, Carlon et al. 
estimate two of four independent critical exponents. De- 
pending on the diffusion rate d, their estimates for = z 
vary in the range 1.60(5) . . . 1.87(3) while (3/v± is found 
to be close to 0.5. Since these values are close the the 
PC exponents z = 1.749(5) and 0/v± = 0.499(2), they 
suggest that the transition might belong to the PC uni- 
versality class. 

The conjectured PC transition poses a puzzle. In all 
cases investigated so far, the PC class requires an ex- 
act symmetry on the level of microscopic rules. In 1+1 
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FIG. 1. Typical space-time trajectories of four differ- 
ent critical nonequilibrium processes starting with a fully oc- 
cupied lattice. The four panels show directed percolation 
(DP), a branching annihilating random walk with two off- 
spring (BAW2), and the pair contact process without (PCP) 
and with diffusion (PCPD). 

dimensions this symmetry may be realized either as a 
parity conservation law or as an explicit Z2 symmetry 
relating two absorbing states. In the PCPD, however, 
the dynamic rules are neither parity conserving nor in- 
variant under an obvious symmetry transformation. Yet 
how can the critical properties of the transition change 
without introducing or breaking a symmetry? As a pos- 
sible way out, there could be a hidden symmetry in the 
model, but we have good reasons to believe that there 
is no such hidden symmetry or conservation law in the 
PCPD. This would imply that the PC class is not char- 
acterized by a "hard" Zi symmetry on the microscopic 
level, rather, it may be sufficient to have a "soft" equiv- 
alence of two different absorbing states in the sense that 
they are reached by the dynamics with the same proba- 
bility. 

In this paper I suggest that the transition in the PCPD 
might belong to a different yet unknown universality 
class. The reasoning is based on the conservative point 
of view that a "soft" equivalence between two absorb- 
ing states is not sufficient to obtain PC critical behavior. 
As described in Ref. [Q, the essence of the PC class is a 
competition between two types of absorbing domain that 
are related by an exact Z2 symmetry. Close to criticality 
these growing domains are separated by localized regions 
of activity. In I +1 dimensions, these active regions may 
be interpreted as kinks between oppositely oriented do- 
mains, which, by their very nature, perform an unbiased 
parity-conserving branching-annihilating random walk. 
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FIG. 2. The density of particles p\(t) and the density 
of pairs p2(t) times i 0,25 as a function of time measured in 
Monte Carlo steps. Upward (downward) curvature indicates 
deviations from the critical point. The dashed line represents 
the slope corresponding to the PC exponent S ~ 0.285. 

In the PCPD, however, it is impossible to give an exact 
definition of "absorbing domains." We can, of course, 
consider empty intervals without particles as absorbing 
domains. Yet, what is the meaning of a domain with 
only one diffusing particle? And even if such a definition 
were meaningful, what would be the boundary between 
an empty and a "single-particle" domain? Moreover, in 
PC models there are two separate sectors of the dynamics 
(namely, with even and odd particle numbers), whereas 
there are no such sectors in the PCPD. In fact, even when 
looking at typical space-time trajectories, the PCPD dif- 
fers significantly from a standard branching-annihilating 
random walk with two offspring (see Fig. |l|). In par- 
ticular, offspring production in the PCPD occurs spon- 
taneously in the bulk when two diffusing particles meet, 
whereas a branching-annihilating random walk generates 
offspring all along the particle trajectories. Therefore, it 
is reasonable to expect that the two critical phenomena 
are not fully equivalent. 

In order to investigate this question in more detail, it is 
useful to compare the DMRG estimates with numerical 
results obtained by Monte Carlo simulations. It is impor- 
tant to note that there are two possible order parameters, 
namely, the particle density 

i 

and the density of pairs of particles 
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where L is the system size and Si(t) — 0, 1 denotes the 
state of site i at time t. Performing high-precision simu- 
lations it turns out that the critical behavior at the tran- 
sition is characterized by unusually strong corrections to 
scaling JpO[ . These scaling corrections are demonstrated 
in Fig. ypwhere the temporal decay of the two order pa- 
rameters for d — 0.1 is shown as a function of time run- 
ning up to almost 10 6 time steps. The pronounced cur- 
vature of the data in the double-logarithmic plot demon- 
strates the presence of strong corrections to scaling. In- 
terestingly, the two curves bend in opposite directions 
and tend toward the same slope. Thus, in contrast to the 
mean- held prediction, pi(t) and pi{t) seem to scale with 
the same exponent. Discriminating between the negative 
curvature of pi(t) and the positive curvature of P2(t), we 
estimate the critical point and the exponent 5 = fijvn by 



p c = 0.1112(1), 5 = ^||=0.25(2). 



(3) 



While this estimate deviates only slightly from the 
known PC value 0.286(2), other exponents deviate more 
significantly. Performing dynamic simulations starting 
with a single pair of particles we measure the sur- 
vival probability P(t) that the system has not yet reached 
one of the two absorbing states p3fl , the average number 
of particles Ni(t) and pairs N 2 {t), and the mean square 
spreading from the origin R 2 (t) averaged over the sur- 
viving runs. At criticality, these quantities should obey 
asymptotic power laws, P(t) ~ t~ s ' , Ni(t) ~ N 2 (t) ~ tf 1 , 
and R 2 (t) ~ t 2 / z , with certain dynamical exponents 8' 
und 77. Notice that in non-DP spreading processes the 
two exponents 8 = (3/v\\ and 8' — [3' may be different. 
Going up to 2 x 10 5 time steps we obtain the estimates 



5' = 0.13(2), 77 = 0.13(3), 2=1.83(5). 



(4) 



Although the precision of these simulations is only mod- 
erate, the estimates differ significantly from the PC ex- 
ponents 8' = 0.286, r] = in the even sector and 8' = 0, 
77 = 0.285 in the odd sector. The exponent z, on the 
other hand, seems to be close to the PC value 1.75. 

The most striking deviation is observed in the expo- 
nent /3, which is not accessible in DMRG studies. Here 
the estimates seem to decrease with increasing numerical 
effort. As an upper bound we find 



(3 < 0.67. 



(5) 



Even more recently, Odor studied a slightly different ver- 
sion of the PCPD on a parallel computer, reporting the 
estimate {3 — 0.58(1) [^2) which is incompatible with the 
PC exponent (3 = 0.92(2). 

In summary the critical behavior of the PCPD is af- 
fected by strong corrections to scaling, wherefore it is 
extremely difficult to estimate the critical exponents. Al- 
though DMRG estimates presented in |17| are very ac- 
curate, they have to be taken with care since they are 



affected by scaling corrections as well. Thus, the appar- 
ent coincidence with the exponents of the PC class may 
be accidental. Comparing other exponents, in particular 
the density exponent (3 and the cluster exponents 8' and 
77, the PC hypothesis can be ruled out. 
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